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Abstract 

Surface detector arrays are designed to measure the spectrum and composition of high-energy cosmic rays by detecting 
the secondary particle flux of the Extensive Air Showers (EAS) induced by the primary cosmic rays. Electromagnetic 
particles and muons constitute the dominant contribution to the ground detector signals. In this paper, we show 
that the ground signal deposit of an EAS can be described in terms of only very few parameters: the primary energy 
E, the zenith angle 8, the distance of the shower maximum Xmax to the ground, and a muon flux normalization 
. This set of physical parameters is sufficient to predict the average particle fluxes at ground level to around 10% 
accuracy. We show that this is valid for hadronic air showers, using the two standard hadronic interaction models 
used in cosmic ray physics, QGSJetll and Sibyll, and for hadronic primaries from protons to iron. Based on this 
model, a new approach to calibrating the energy scale of ground array experiments is developed, which factors out 
the model dependence inherent in such calibrations up to now. Additionally, the method yields a measurement of 
the average number of muons in EAS. The measured distribution of A*'^ of cosmic ray air showers can then be 
analysed, in conjunction with measurements of Xmax from fluorescence detectors, to put constraints on the cosmic 
ray composition and hadronic interaction models. 
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1. Introduction 

The origin of Ultra-High Energy Cosmic Rays (UHECR, with energies E > 10^* eV) still remains a 
mystery. Experimental results [1-3] suggest that the UHECR flux is composed predominantly of hadronic 
primary particles. As charged particles, they suffer deflections in cosmic magnetic fields and do not point 
back directly to their sources. Indirect proofs of their origin are necessary instead: the precise measurement 
of the energy spectrum, an estimation of the mass composition and its evolution with energy, and angular 
anisotropies are the three main handles on disentangling this almost ccntury-old problem. Due to the inter- 
action with the cosmic microwave background, UHECR suffer energy losses which limit their propagation 
distance [4,5]. This "GZK horizon", indications of which have already been observed in the UHECR spec- 
trum [6-8] , depends very sensitively on the energy and mass of the cosmic ray. In order to discern between 
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different source scenarios, and to disentangle source characteristics from the effects of propagation, a precise 
knowledge of the energies of UHECR is crucial. Constraints on the composition of the cosmic ray flux at 
the highest energies will supply additional fundamental insight. 

Due to the low fluxes at ultra-high energies, the detection of UHECR can only be achieved by mea- 
suring Extensive Air Showers (EAS), cascades of secondary particles resulting from the interaction of the 
primary cosmic rays with the Earth's atmosphere. The measurement of the cosmic ray energy, flux, and 
mass composition relies on an understanding of this phenomenon. 

Two main EAS detection techniques were developed over the years (see [9] for a review): surface detectors 
(SD) detect the particle flux of an EAS at a particular stage of the shower development; fluorescence detectors 
(FD) measure the shower development through nitrogen fluorescence emission induced by the electrons in the 
shower. The modeling of EAS through Monte Carlo simulations is needed in both fluorescence and surface 
detector experiments in order to interpret the data. We will show that hadronic EAS can be characterized, 
to a remarkable degree of precision, by only three parameters: the primary energy the depth of shower 
maximum A"niax, and an overall normalization of the muon component, which we call N^. This is what we 
will call air shower universality [10]. The parameters Amax and are linked to the mass of the primary 
particle, ranging from proton to iron, and are subject to showcr-to-shower fluctuations; proton showers have 
a larger depth of shower maximum than iron showers, while iron showers contain ~ 40% more muons than 
those induced by protons. Once measured, and Amax have to be compared with simulations to infer 
the cosmic ray composition and place constraints on hadronic interaction models. Previous studies have 
demonstrated that the energy spectra and angular distributions of electromagnetic particles [11,12], as well 
as the lateral distribution of energy deposit close to the shower core [13] are all universal, i.e. they are 
functions of E, Amax, and the atmospheric depth A only. ^ For studies of shower universality in the context 
of ground detectors, see [14-16]. EAS induced by photons show somewhat different properties, due to the 
absent hadronic cascade. Hence, it remains to be investigated to what extent the hadronic EAS universality 
studied here applies to photon showers. 

By sampling the longitudinal development of the electromagnetic shower component close to the core, 
fluorescence detectors measure both Amax and E. The systematic uncertainty in the energy E is typically 
25%, mainly due to the uncertainties in the air fluorescence yield. A surface detector only samples the 
properties of an EAS at a given stage of the shower development and at several points at different distances 
r from the shower axis. Rather than using the signal integrated over all distances, a quantity which shows 
large fluctuations, Hillas [17] proposed to use the signal at a given distance r from the shower axis, S{r), 
as a measure of the shower size, connected with the primary energy. The distance where experimental 
uncertainties in the size determination are minimized (the optimal distance ropt [18]) is mainly determined by 
the experiment geometry, i.e. the spacing between surface detectors. S'(ropt) is then related with the primary 
energy of the incoming cosmic ray using Monte Carlo simulations. This calibration has large systematics 
due to uncertainties in the hadronic models and the unknown primary cosmic ray composition. 

In this paper, we will show how to use air shower universality to determine the calibration of a surface 
detector in a model-independent way. The signal <S'(ropt) is the sum of two components: an electromagnetic 
part which is well-understood and to a good approximation depends only on E and Xmax of the shower; and 
a muon part which, in addition to E and Amax, depends on the model and primary composition in terms 
of an overall normalization. The muon fraction can be determined by requiring that the shape of the zenith 
angle dependence of 5(ropt) at a fixed energy, which depends on the muon normalization A^j^, match the 
observed one. This method determines the energy scale of the experiment as well as the average number of 
muons produced in the air showers at a given energy. 

Subsequently, we will apply air shower universality to data collected by a hybrid experiment, which com- 
bines the fluorescence technique with a surface detector. In this case, the calibration of the surface detector 
can be done almost independently of hadronic models and composition by using a small subset of the data 
(hybrid events) which are simultaneously measured by the fluorescence and the surface detector. Applying 
our method to hybrid data yields an event-by-event measurement of the muon content of the shower. This 
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can be used as an independent cross-check of the measurement from the surface detector alone. Since the 

electromagnetic contribution to the signal varies with zenith angle, a hybrid measurement of A'^; at different 
zenith angles probes whether the electromagnetic part is described correctly by simulations, a key ingredi- 
ent in our study. Conversely, the surface detector energy scale obtained with the universality-based method 
offers a cross-check of the hybrid calibration of the surface detector, which uses the fluorescence energy 
measurement. 

In this work, we will not use data from any experiment, but wc will use the Pierre Auger Observatory as 
a case of study. First results from this method applied to Auger data have already been presented in [19]. 
While we adopt the specifications of this experiment, the method presented here can be applied to any other 
surface detector (for example, AGASA [20]) or hybrid experiment (for example. Telescope Array [21]). The 
paper is organized as follows: in Sec. 2 we explain air shower universality, verifying it with the two standard 
high energy hadronic models used in cosmic ray physics (QGS.JetII and Sibyll); in Sec. 3 the limits of air 
shower universality are shown; Sec. 4 presents the; nicrthod of obtaining from the surface detector and 
determining the surface detector energy scale: in Sec. -5 wc validate the method using a simple Monte Carlo 
approach; Sec. 6 shows how the approach can be applied to hybrid events; finally, the application to other 
experiments is discussed in Sec. 7; we conclude in Sec. 8. 

2. Extensive air shower universality at lairge core distances 

The results presented in this paper were obtained from a library of simulated EAS. We used CORSIKA 
6.500 [22] with hadronic interaction models QGSJetll [23,24] and Fluka [25] (proton and iron primaries, 
energies lO^^-^ - 10^" eV), and Sibyll [26,27] / Fluka as well as QGSJctll / Gheisha [28] (proton at lO^^ eV). 
For each primary/energy combination, we simulated 80 showers each at 7 zenith angles ranging from 0° to 
60°. Statistical thinning was employed in the simulations as describeds in [29], at a thinning threshold of 

£ = 10-6. 

Using lookup tables generated with GEANT4 simulations [30], we calculated the average response of a 
cylindrical water Cherenkov detector (height 1.2m, cross section 10 m^, similar to the type used in the Pierre 
Auger observatory) to each shower particle hitting the ground. See Sec. 7 for a discussion of the applicability 
to other experiments. The signals were calculated in two difi^erent approaches: 1.) Ground plane signals: The 
response is calculated for a realistic water tank on the ground. 2.) Shower plane signals: The response is 
calculated for a fiducial flat detector (with the same average particle response as the water tank) placed in 
the plane orthogonal to the shower axis {shower plane). Signals calculated in the shower plane procedure are 
not affected by detector geometrical effects, and therefore independent of the zenith angle. For details on the 
signal calculation, see appendix A. The shower plane signals will be useful to verify air shower universality, 
while the ground plane signals will be needed for the application to a realistic experiment (in our case the 
Pierre Auger Observatory). 

Due to the statistical thinning procedure employed in the shower simulation, particles were collected in a 
sampling area of width 0.1 in log^Q r centered around the shower core distance r considered. This ensures, 
for a wide range of slopes of the lateral distribution, that the median radius of the energy deposit in the 
sampling area is indeed r. Signals are calculated in 18 azimuthal sectors, and normalized relative to the 
signal deposited by a vertically incident muon (VEM), a standard practice in surface detectors using the 
water Cherenkov technique. 

To describe the stage of the shower development, we use the variable DX, defined as the distance from 
the detector to the shower maximum measured along the shower axis (in g/cm^). For a tank on ground at 
a distance r from the shower axis, DX is: 

DX = Xgr sec 6 — Xjaax — f COS C tan 6 Pair (1) 

where Xgr is the vertical depth of the atmosphere, is the zenith angle, X^ax is the slant depth of shower 
maximum, and C, is the azimuthal angle in the shower plane such that C = corresponds to a tank below 
the shower axis, pair ~ 10~^g/cm^ is the density of air at ground level (see also Fig. A.l in the appendix). 
Often, we will consider signals averaged over azimuth. In this case, DX is simply given by: 
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Fig. 1. Left panel: Simulated electromagnetic shower plane signals at r = 1000 m for proton (red dots) and iron showers (blue 
circles) at 10^^ eV as a function of DX. The showers are simulated with QGSJetll/Fluka at discrete zenith angles spanning 
0° to 60°. Right panel: Simulated electromagnetic shower plane signals vs. DX for different primaries and hadronic models, 
relative to the prediction for proton showers when using QGSJetll/Fluka. 



DX = X„ sec e- Xr. 
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2.1. Electromagnetic and muon shower plane signals 



At large core distances (r > 100 m), the particle flux of EAS at ground is dominated by electromagnetic 
particles (e+, e^, 7) and muons. Throughout the paper, we include the signal from the electromagnetic 
products of in-flight muon decay in the muon contribution, separating it from the 'pure' electromagnetic 
part. 

Fig. 1 (left panel) shows the electromagnetic shower plane signals Sem of simulated proton and iron showers 
at 10^^ eV as a function of DX (in g/cm^) for a core distance of 1000 m. For each tank, we calculate the 
corresponding DX via Eq. (1) using the azimuth angle of the tank. Since zenith angle dependent detector 
geometry effects are removed in the shower plane treatment, we are able to compare the signals from a wide 
range of zenith angles. The electromagnetic signal shows a strong evolution with DX, reaching a maximum 
at DXpcak and rapidly attenuated for larger DX. DXpcak depends on core distance, being Og/cm^ very close 
to the core and w 200g/cm^ at 1000 m. This shift is only mildly dependent on r in the range 400— 1600 m 
and it can be naturally explained by diffusion of electromagnetic particles away from the shower axis. 
Note that the overall electromagnetic signal as well as its evolution arc slightly different for protons and 
iron. This is apparent in the right panel of Fig. 1, where the ratio of the signals obtained from different 
primary/model combinations to proton-QGSJetll is shown as a function of DX. The differences between 
models are around 5-10%, smaller than the deviation between proton and iron. This result is an extension 
to large r of previous results [13,12,31] on the universality of the electromagnetic EAS component at small 
core distances. We address the difference (~ 15%) in Sem between protons and iron in Sec. 3. 

Fig. 2 (right panel) shows the evolution of the muon signal S*^; with DX (again for proton and iron showers 
at i? = 10^^ eV and r=1000 m). shows a distinctly different behavior: it peaks at DXpcak ~ 400 g/cm^, 
and it is attenuated much more slowly than 5em- As expected, there is a dependence of the absolute 
normalization of the signals on the primary particle and hadronic model, which is clearly seen in Fig. 2 
(right panel) where we again show 5*^ obtained for different primaries and models relative to that of proton- 
QGSJetll. As for S'em, only differences in normalization and not in shape are apparent. We verified that the 
primary- and model-indcpendcncc of the electromagnetic and muon signal evolution holds for shower core 
distances between 100 m and 1000 m. 
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Fig. 2. Le/t panel: Simulated muon signals at r = 1000 m vs. DX for the same 10^^ eV proton (red dots) and iron showers 
(blue circles) as in Fig. 1. Note that Sfj, includes the contribution from muon decay products. Right panel: Simulated muon 
signals vs. DX for different primaries and hadronic models relative to the muon signal predicted for proton showers when using 
QGSJetll/Fluka. Note the difference in scale compared to Fig. 1. 



2.2. Ground signal parameterization 

In the previous section, we have shown that the evolution of the shower plane signals at a given shower 
core distance is only very weakly dependent on the primary particle or hadronic model considered. Therefore, 
a simple parameterization of the signals is possible. In this work, since we use the Pierre Auger Observatory 
as a case of study, we perform such a parameterization for r=1000 m. It has been shown [32] that the main 
observable in the surface detector of the Pierre Auger Observatory, S'(IOOO), is indeed a good measurement of 
the azimuth-averaged signal of particles at a core distance of r = 1000 m. Hence, we separately parameterize 
the azimuth-averaged ground plane electromagnetic and muon signal at 10^^ eV (5'em(1000) and 5^(1000), 
the total predicted signal being the sum of both), using the incomplete gamma, or Gaisser-Hillas-type 
function: 

5(1000, DX) = 5_ ^-^^^-—^^ j exp j , a = (3) 

The four free parameters of this function are: S'max (the peak signal at 1000 m); Z?Xpcak(the slant depth 
relative to the overall shower maximum where the peak signal is reached); A (the attenuation length after 
the maximum); and Xq (an additional shape parameter). 

Fig. 3 shows the results of the fit for the muon signal. We have simultaneously fitted the predictions for 
different primaries (proton, iron) and different models (QGSJetll, Sibyll), keeping a separate normaliza- 
tion (5max) for each, while A and DX-p^ak are common to all. Xq is not fitted but fixed to —200 g/cw? . 
The resulting parameters are summarized in Table 1, with 5max given for proton-QGSJetll. For the other 
model/primary combinations, we define a relative muon normalization given by iV^ = >S'max/'5'max;rcf , where 
we take proton-QGSJetll as the reference S'maxirof . The A^^ for different models and primaries are listed in 
Table 2. 

In the case of the electromagnetic signal, we have to take into account the detector geometrical effects, 
which cause differences in the signals from showers at two different zenith angles with the same DX (see 
appendix A). Hence, we have to find a parameterization for S-em{DX^ 9). The first step is to parameterize, 
for each of the 7 simulated zenith angles, the dependence of 5em on DX; a linear function is found to be 
sufficient due to the limited DX range at a fixed 9. We scaled the proton and iron signals by 1 -I- a and 1 — a, 
respectively (with a < 0.06), to account for the deviations from universality. The deviations in the ground 
plane signals arc slightly smaller than those shown in the previous section for the shower plane signals. Fig. 4 
shows the results of the fits together with the direct Monte Carlo results, for the 7 fixed values of zenith 
angle. 
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Fig. 3. Parameterization of the muon ground plane sig- Fig. 4. Electromagnetic ground plane signals (proton— dots; 
nal (E = IQl^ eV) at r = 1000 m using Gais- iron— crosses; QGSJetll) in zenith angle bands (color-coded). 
ser-Hillas functions (Eq. (3)). Red dots (crosses) denote pro- Proton and iron signals have been scaled symmetrically. A 
ton-QGSJetll (proton-Sibyll) showers, blue circles (asterisks) linear function is fit to the signal separately for each zenith 
are iron-QGSJctll (iron-Sibyll). The normalization is left free angle band, 
for each model/primary combination, while the other param- 
eters are common to all. 

In the second step, we fit a Gaisser-Hillas-type function (of DX) to Sem{DX, 9) for each X^ax considered. 
The Gaisser-Hillas function is fitted to 7 equal-weight data points "predicted" from the 7 Hnear fits of the 
first step. This is equivalent to a parameterization of the dependence of ^em on at a fixed X„iax, using an 
intermediate variable DX{Xa^s.^,9). Table 1 gives the results for X^ax — 750g/cm^ and Xgr = 875g/cm^. 
The second step may be applied to any value of Xmax, yielding a continuous function 5'em(^, ^-^(-^max, d)) 
which depends on 6 both explicitly and implicitly via DX . By contrast, since muons deposit a signal which is 
proportional to their pathlength in the water tank, the tank acts as a volume detector: the smaller projected 
area at higher zenith angle is canceled by the longer average tracklength. Hence, the average muon signal 
Sfi{DX) does not show an explicit 6 dependence. 

At a fixed energy (here, 10 EeV), the parameterization presented above determines the average ground 
signal of a shower (at r = 1000 m, azimuth-averaged): 

5(1000) = 5em(0, i^^(^max, e)) + N,, ■ S^,.rc{{DX{X^,^, 9)) (4) 

Here, Sem denotes the parameterized electromagnetic signal, and S^-ref is the reference muon signal which 
we take to be proton-QGSJetll. Hence, there are only three free parameters describing the average shower 
at this energy: the zenith angle 9; the depth of shower maximum Xmax! and the normalization of the muon 
signal Nfj, (relative to proton-QGSJetH). 

We used the library of proton and iron showers with energies of 10^^ — 10^" eV to investigate the energy 
dependence of the evolution of S'em and 5^ with DX. The electromagnetic signal normalization shows 
an energy scaling of 5'max;EM cx; E^-^'^ (see also Sec. 3), while for the muon signal S'niax;/^ oc with 
a = 0.9 . . . 0.95, depending on the hadronic model. All other fit parameters in Eq. (3) are independent of 
the primary energy in this energy range, for both S'em and S^, to within 5%. 

Hence, Eq. (4) can be straightforwardly extended to other energies: 
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Table 1 

Fit parameters of the Gaisser-Hillas parameterization (Eq. (3)) of the universal electromagnetic and muon sig nal at 10^^ eV. 
The electromagnetic parameterization is for a fixed Xmax and Xgr = 875 g/cm^. For 5^ max, the value of proton-QGSJetll is 
given (for the other primaries and models relative to proton-QGSJetll, see Table 2). 
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Fig. 5. Left panel: Distribution of the relative deviations of the electromagnetic ground plane signals of showers at 10^^ eV 
(QGSJetll) from the parameterization (Sec. 2.2). The red (solid) lino is for proton, while the blue (dashed) is for iron. Right 
panel: The same for the muon ground plane signals. 



S(1000, E) = 5em(10 EeV, 6, DX{X,^,^, 6)) 
+ N^{E) ■ 5^;„f (10 EeV, DX{X 
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As the energy scaling of 5^ is slightly model-dependent, we treat it as an unknown and define Nf^{E) as the 
muon normalization at the energy E with respect to the proton-QGSJetll reference at the fixed energy of 
10 EeV. 



2.3. Shower fluctuations 

In addition to the overall behavior of the signals with DX parameterized before, both electromagnetic 
and muon signals show fluctuations around the mean value. Fig. 5 shows the relative deviations of the 
ground plane 5em (left panel) and (right panel) from the parameterization for proton and iron showers 
(10^^ eV, QGSJetll). These distributions contain showers from all zenith angles; no dependence of the 
relative fluctuations on zenith angle has been found. Note that the proton and iron electromagnetic signals 
are slightly shifted from due to the universality violation (Fig. 1), whereas the deviations are centered 
around for the muon signals (we used the corresponding muon signal normalizations for proton/iron). 

The spread of the distribution shown in Fig. 1 and Fig. 2 has a contribution from the artificial fluctuations 
due to the thinning procedure used in the simulations. The fluctuations due to thinning can be estimated 
from vertical showers: since we expect the same signal in all azimuth sectors, thinning fluctuations are 
expected to be the dominant source of the variance between sectors. We find (Xthin — 6.5% for S'em and 4% 
for Sfj_. We can then subtract the uncorrelated thinning variance from the total signal fiuctuations to obtain 
the shower-to-showcr fiuctuations. Note that since we compare shower signals with the parameterization of 
the average signal at the same distance to ground DX, the signal fiuctuations shown here are not caused 
by the fiuctuations in the depth of shower maximum. The latter ones will induce additional fluctuations 
(mainly in the electromagnetic signal) that can be straightforwardly calculated convolving the fluctuations 
in Xniax with the signal parameterization. 

We also parameterized the distribution of Xmax for different primaries and models, using the following 
functional form: 

— cx X e , < a; < oo; x = ^ ^ + 5, (6) 

"A max Tx 

where (Xmax) denotes the mean depth of shower maximum, and tx is related to the RMS of the distribution 
via Tx — RMS(Xniax)/\/5- This asymmetric distribution is found to be a good fit to the Xmax distributions 
for different primaries and models. 
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Table 2 summarizes the magnitude of fluctuations in Sem and as well as tx for protons and iron 
using different hadronic models. Clearly, the shower-to-shower fluctuations in signal as well as Xmax are 
independent of the hadronic model considered, but depend quite strongly on the primary particle. Hence, if 
measured, fluctuations can serve as a robust, model-independent indicator of composition. Our simulations 
predict that these fluctuations depend only very weakly on energy. 



3. Limits of universality 

The main discernible deviation from the universality approach adopted here is the difference in electro- 
magnetic signal between proton and iron showers. This difference, which we refer to as universality violation, 
is larger than the differences found in the overall energy deposit in the atmosphere for proton and iron show- 
ers (for which in fact one finds the opposite effect: the so-called missing energy is larger for iron showers, 
[33]). Since we include muon decay products in the muon signal, this deviation is unrelated to the differences 
in muon content between protons and iron. 

Fig. 6 shows the ratio of the number flux of electromagnetic particles for different combinations of pri- 
maries/models to the reference (proton/QGSJetll) as a fmiction of DX (again r=1000 m and E = 10^^ eV). 
The differences between protons and iron is much smaller than in the case of the signals (Fig. 1), pointing 
to a slightly harder energy spectrum for electromagnetic particles at large r in iron showers compared to 
proton showers. We also found that the discrepancy becomes smaller at smaller core distances. We have 
verified that the differences are independent of the details of nuclear fragmentation of the primary iron 
nuclei. This means that the nuclear binding of the 56 nucleons is not important for the EAS development. 
In other words, the superposition model holds, i.e., an iron shower can be considered as a superposition of 
56 proton showers at 1/ 56th of the primary energy. 

This implies that the universality violation is due to a violation of strict linear energy scaling of the 
electromagnetic signal in hadronic shower simulations: if the electromagnetic signal scales as E", a < 1, 
then the signal of an iron shower will be a factor of 56^~" times larger than that of a proton shower at the 
same energy. In order to explain the observed difference of ~ 15% in the shower plane signals, we would 
infer a ~ 0.97. We have parameterized the electromagnetic signals for different energies and indeed foimd 
that the amplitude ^max of the signal (Eq. (3)) scales as E^-^"^ at r = 1000 m, with a approaching 1 as 
r — » (Fig. 7). Note that, by parameterizing the complete evolution of the signal with DX, we take out the 
effects of the energy dependence of Xmax- 

This violation of perfect energy scaling of the electromagnetic signal can be due to several reasons. The 
injection rate of energy into the electromagnetic part via 7r° decay as well as the energy spectrum of secondary 
71-0 might evolve with primary energy. In additicm, the NKG theory of pure electromagnetic showers also 
predicts a slight deviation from perfect energy scaling of the particle flux on ground. These effects are 
currently under investigation. 





RMS(5bm)/Sem 


RMS(5^)/5^ 


(Xmax> (10 EeV) 


Tx 




Proton 












QGSJetll 


7.9% 


11.8% 


787.8 g/cm^ 


25 A g/cm^ 


1 


Sibyll 


8.9% 


12.4% 


795.8 g/cm2 


24.1 g/cm^ 


0.87 


Iron 












QGSJetll 


5.4% 


3.5% 


708.7 g/cm^ 


10.9 g/cm^ 


1.40 


Sibyll 


4.8% 


4.0% 


696.5 g/cm^ 


10.2 g/cm^ 


1.27 



Table 2 

Relative shower-to-shower fluctuations of the electromagnetic and muon signals and parameters of the Xmax distribution 
derived from QGSJetll and Sibyll showers at lO^^ eV. The muon signal normalization A^^ relative to proton-QGSJetll for 
the different models is also shown. Note the differences in the absolute value of (Xmax) and Nfj,, while the fluctuations are 
model-independent. 
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Fig. 6. Number flux of electromagnetic particles in the shower 
plane at r = 1000 m for different primaries and hadronic 
models at 10^^ eV, relative to that of proton-QGSJetll. 



Fig. 7. The parameter Smax of the parameterization Eq. (3) of 
electromagnetic shower plane signals as a function of primary 
energy, for different core distances r. Power-law fits and the 
resulting exponents are indicated. 



4. Determining the muon normalization using the constant intensity method 



One of the main challenges of a cosmic ray surface detector is to convert the ground signal S{r) to a 
primary energy. As mentioned in Sec. 2.2, the universality-based signal parameterization has three free 
parameters. Apart from the zenith angle 9 which is well measured along with the signal [32], the depth 
of shower maximum ATmax and muon normalization N^^ (with respect to the reference signal at the fixed 
energy of 10 EeV) remain to be determined. Once these are known, Eq. (5) provides a one-to-one mapping 
of ground signal and energy, i.e. a model-independent energy scale of the experiment. 

The mean depth (Amax) of showers has been measured as a function of energy from experiments using 
the air fluorescence technique, e.g. HiRes [34] and Auger [35]. The knowledge of (Amax) is important as it 
determines the average distance to shower maximum DX for a given zenith angle, and the electromagnetic 
signal evolves strongly with DX. The overall precision of these (Amax) measurements is better than 20g/cm^, 
and this small uncertainty in (Amax) has only a limited effect on the estimated electromagnetic signal. Fig. 8 
shows the limited effect of varying Amax by ±14g/cm^ (the current measurement uncertainty at 10 EeV 
reported by Auger [35]) on Sem as a function of sec6'. 

The main uncertainty in determining the energy scale of surface detectors is thus in the value of N^. 
Fortunately, one can make use of the different behavior of Sem and with DX (and hence, sec 9) to 
measure via the constant intensity method (Fig. 9): dividing the data set into equal exposure bins in 
zenith angle, i.e., bins of sin^ 9, a correct signal-to-energy convertor should yield the same number of events 
in each bin with measured signal greater than the parameterized signal at a fixed energy. This is due to the 
isotropy (^-independence) of the cosmic ray flux, which requires that the number of events N{> E) above 
a fixed energy E should be equal in equal exposure bins. 

Fig. 9 (upper panel) shows the zenith angle dependence of the signal (Eq. (4)) for a fixed energy of 10^^ 
eV and different values of A"^. Apart from the overall change in signal, it is evident that the smaller the 
A"^, the steeper the 9 dependence is. We now divide a simulated ground detector data set with a "true" 
Af^(lO^^eV) = 1 (see Sec. 5 for details on the simulation) in equal exposure bins in zenith angle. Given 
a muon normalization, we calculate the number of events in each bin that are above a given reference 
energy (here i?rof=10^® eV), according to Eq. (4) with the given A'^. We then adjust Ni^i{Ei-ci) in the signal 
parameterization Eq. (4) to the value which gives an equal number of events A^(> S^Ej-d, 9)) in each zenith 
angle bin (lower panel in Fig. 9). Clearly, a too low value of A"^ results in an excess of events at high 9 (the 
parameterized signal has a too steep attenuation with sec 9), whereas a too high A^^ results in a deficit of 
high zenith angle events {sec 9 attenuation too shallow). In this calculation we used the (Amax) at 10 EeV 
reported in [35] to calculate the ground plane signals. Note that an A"^ of 1.1 gives a flat distribution, 
whereas the "true" A"^ used is 1.0. This bias in the A'^ measurement will be addressed in the next section. 
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Fig. 8. Parameterized electromagnetic signal at 10^ eV vs. 
sec 6 with Xmax = 750g/cm^ (black line, the dots indicate the 
signal parameterized at each zenith angle). The shaded band 
shows the effect on Sbm of a variation of Jfmax by ±14g/cm^. 



Fig. 9. Upper panel: the signal parameterization at 10 EeV 
Eq. (4) vs. sec d for different Nf^ (black/solid— 1.1, 
red/dashed— 0.5, blue/dash-dotted— 2.0). Lower panel: his- 
tograms of number of events above the parameterized signal 
in equal exposure bins, obtained from a Monte Carlo data set 
with a true A^^ of 1, for the different A^^ values shown in the 
upper panel. 



For a range of values, we then calculate the x^/dof of the event histogram relative to a flat distribution 
in sin^ 9. Fitting a parabola to the function x^i^fj.) yields the best-fit iV^et and its error (Tat^: 

xHN,)^xL.+ (^^^^^y (7) 

For a data set comparable to current Auger statistics (~11 000 events above 3 EeV [7]), we expect 
a statistical error of ctn^ = 0.1. Once iV^ is known, the knowledge of S'em (within the uncertainty of 

±6% due to universality violation) determines a model-independent energy scale, with a statistical error 
of ctaTj, ■ <S'^;rcf around 4%. The constant intensity method can be extended to other energies, using the 
energy-dependent parameterization Eq. (5) in Sec. 2.2. This yields a measurement of N^{E), comparable 
to the measurement of (X,„ax) in its sensitivity to the primary composition. Table 3 contains a summary of 
the expected statistical and systematic errors from current and upcoming experiments. 

In Fig. 10 we show possible results of this measurement, the integral measurement of Nfj_{E) (solid black 
line, corrected for the bias, see Sec. 5) and with la statistical error band (shaded) after a three year Auger 
exposure. Here, we took Nfj_{E) = 1.2(£'/10 EeV)"'^^ as fiducial value. As the cosmic ray spectrum drops 
rapidly with energy (oc E~^), the average energy of cosmic rays above a given energy threshold is very close to 
that threshold. Hence, for slow changes of the cosmic ray composition with energy, the Nf^ value determined 
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Fig. 10. The measured muon normalization Nfj, as a function of energy (thick black line) with statitistical error band expected 
from a three-year Auger exposure (shaded), for a fiducial Nfj,{E) = 1.2 (E/10 EeV)°-**5. Also shown are model predictions 
for iron (upper two lines, blue; solid— QGSJetll, dashed— Sibyll) and proton showers (lower two lines, red; solid— QGSJetll, 
dashed-Sibyll). 

from the constant intensity method wiU reflect the actual average value of for cosmic rays at that energy 
(this will be shown in the next section). In case of an abruptly changing composition, the measured N^(E) 
will clearly show evidence for this. However, the interpretation of the integral Nf^ measurement in terms 
of composition will have to rely on a modeling of the composition evolution in this case. In addition, the 
measurement of Nf^{E) can place constraints on hadronic models, whose predictions are shown as lines in 
Fig. 10. 

5. Validating the constant intensity method 

To benchmark and validate the determination of via the constant intensity method, we simulate 
realistic data sets based on our parameterization of the ground signal (Sec. 2.2) and its fluctuations (Sec. 2.3). 
The fluctuations in signal as well as X^ax could have an impact on the measurement of N^^ , since only the 
average values are used to infer N^^ (Sec. 4). Additionally, a mixed composition of the cosmic ray beam could 
bias the measurements. The purpose of this section is to quantify systematic uncertainties of the method 
described before. 

The calculation of a simulated data set proceeds as follows. Event energies are drawn from a spectrum 
dN/dE cx E^^-^ in the range 10^^'* — 10^°'^eV, while the zenith angle is drawn from an isotropic distribution 
{0 < 70°). The primary particle type (proton or iron) is chosen at random according to a given mixture. 





Muon normalization N/j, 


Energy scale 


Statistical error 






current Auger 


0.1 


4% 


Systematic errors 






(^max) uncertainty (14g/cm2 [35]) 


+0.05 / -0.07 


-fO.5% / -2% 


Universality violation 


+0.01 / -0.04 


-f 3% / -4% 


bias 


<10% 


<5% 



Table 3 

Expected statistical (for current Auger exposure) and systematic errors on the muon normalization N^j, and energy scale 
S{e = eo, E = 10 EeV) at 10 EeV. 
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Fig. 11. Constraints in the TV^-energy scale plane placed by three independent measurements: the constant intensity method 
(black dot with error); vertical hybrid events (blue solid line with shaded error band); inclined hybrid events (red dashed 
lines). An error of 25% in the fluorescence energy scale is indicated (vertical lines). All values are calculated for a three year 
Auger-equivalent exposure. The fiducial A^^ is 1.2. 



The depth of shower maximum is then drawn from the distribution Eq. (6) with the parameters for the 
given primary (we adopt the parameters from QGSJetll; this has no influence on our conclusions). An N^, is 
determined according to the primary. With E, 9, Xmax, and Nf^ given, the ground signal can be determined 
via Eq. (4) (we scale Sem with and with E^-^). The two signal components are fluctuated according 

to the primary (see Table 2). Finally, we cut events according to a simple trigger depending on the ground 
signal, and apply signal reconstruction uncertainties as reported by the Auger observatory [32]. The main 
characteristics of the reconstruction of S'(IOOO) are that it is unbiased at large signals 5(1000) > 10 VEM, 
and that bias and variance increase quickly for signals below 10 VEM. 

A large set of simulated data sets showed that the error calculation according to Eq. (7) is a good estimator 
for the variance of the measurement. However, we found that the constant intensity method yields a 
systematic shift to higher Nf^^ values of about 5-10%. This can be explained by trigger effects and fluctuations. 
Due to the attenuation of the signal with zenith angle at a flxed energy, the resolution gets worse at large 
zenith angles. Additionally, upward fluctuations above the trigger threshold are more important at high 
zenith angles. These two effects, in the presence of a steep spectrum, produce a zenith angle-dependent 
enhancement of the number of events reconstructed above a given energy. This tends to flatten the constant 
intensity curve (Fig. 9), which leads to a higher estimated A^^ value. The bias in is mainly determined 
by the experimental resolution and trigger effects. It depends slightly on the primary composition, ranging 
from 4% for pure iron to 8% for a mixed composition, due to the differing magnitudes of shower-to-shower 
fluctuations. The unknown composition, and imperfect knowledge of the experimental characteristics, lead 
to an uncertainty in the bias which should be included in the systematic error on N^. We assume that this 
error will be smaller than the absolute value of the bias, hence <10%. 

Further systematics of are the violation of universality in the electromagnetic ground signal {'^ ±6%, 
Sec. 2.2), and the uncertainty in the value of (Amax), which translates into a further uncertainty in the 
electromagnetic signal. Table 3 gives a summary of the systematic errors derived for A^^ and the energy 
scale (i.e., S{6 = Oq, E ^ 10 EeV), we take 6*0 ~ 38°, close to the median of the isotropic cosmic rays), all 
evaluated at 10 EeV. We stress that since the universality violation is smaller closer to the shower core, this 
method exhibits signiflcantly smaller systematics when applied to a surface detector measuring the signal 
at r = 600 m, instead of 1000 m as considered here. 
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6. Cross-checks and independent hybrid measurements 

The constant intensity method described above is independent of the primary composition and hadronic 
interaction models (within systematics) , and also independent from other energy calibrations (e.g., fluores- 
cence telescopes). However, it relics on a good understanding of the electromagnetic part of air showers as 
well as detector bias and resolution (Sec. 5). Hence, it is desirable to cross-check the results of the method 
with independent data. 

Hybrid experiments, which simultaneously measure the ground signal as well as fluorescence energy and 
-''^max on an event-by-event basis, allow several cross-checks of the measurement. Due to the uncertainty 
in the energy scale of the fluorescence detector (~ 25%), we introduce a scaling factor (/fd) of the measured 
fluorescence energy. 

Using the parameterized electromagnetic signal Se,m{DX{9, Xjna,x),0, E) (Sec. 2.2), with E = JfdEfd 
and Xmax given by the fluorescence measurement, we can, for a single event, determine the muon signal S/j, at 
a fixed core distance by subtracting the electromagnetic component from the total signal (see Eq. (5)). The 
muon normalization is then given by = S^/S^-rei, where S^-rei is the reference parameterized muon signal 
(proton-QGSJetll) at 10 EeV. Additionally, we can divide the hybrid data into two sets: vertical and inclined 
events, with zenith angles smaller and larger than 60°, respectively. In inclined events the electromagnetic 
signal at ground is essentially negligible, and the ground signal allows for a direct measurement of the muon 
signal. We can then calculate the mean measured (A^^) for vertical and inclined events as a function of fpD- 
Fig. 11 shows the results for vertical (blue line with shaded error band) and inclined (red lines) events for 
an experiment with 3 years Auger-equivalent exposure. In order to measure (Nfj) with hybrid events we 
clearly need to constrain /fd- The black dot in Fig. 11 corresponds to the result of the constant intensity 
method described in the previous sections. It constrains as well as the energy scale. The fluorescence 
energy scale and its current uncertainty {Jfd = 1-0 ± 0.25) are indicated in the graph by vertical lines. 

The crossing point of the three A^^ measurements is an important cross-check of the universality-based 
method: only a correct description of the evolution of the electromagnetic and muon groimd signals will 
lead to a unique crossing point. The value of Jfd that corresponds to the crossing point is a quite powerful 
fluorescence-independent measurement of the energy scale. The statistical uncertainty of this measurement is 
much smaller than the current uncertainty on the fluorescence energy scale, and thus provides for a sensitive 
cross-check. At the same time, experimental efforts to reduce the systematic fluorescence energy uncertainty 
are in progress [36]. 

Hybrid events offer several further ways to place constraints on hadronic models. For example, since 
Sn and DX are measured independently for each hybrid event, the behavior of Sn{DX) can be inferred, 
which contains information on the energy spectrum of muons in UHE air showers. In addition, the measured 
fluctuations of TV^ allow for model-independent constraints on the primary composition, if the reconstruction 
uncertainties are well understood. In addition, observed anomalous values in single events can be used 
to search for non-hadronic primaries. Photon showers have a muon component of about 1/lOth of a proton 
shower and thus could leave a distinctive signature in the observed N^. However, the sensitivity of this 
method to photon showers remains to be studied quantitatively with simulations of photon showers. 

7. Application to other experiments 

The methodology and results presented so far have been specialized to the case of Auger. We now discuss 
the applicability to other current and future experiments. The main experiment characteristic determining 
this application is the ratio of the electromagnetic and muon contributions to the shower size observable 
(5(1000) in the case of Auger). 

If the muon contribution is very small, the signal becomes essentially independent of N^, and only the 
knowledge of (Xmax) is needed to predict the average signal in a model-independent way. An experiment 
operating in this regime (such as AG AS A [37]) is able to experimentally verify the electromagnetic signal 
predicted by simulations. 

Conversely, if the muon contribution dominates even at small zenith angles, the attenuation (i.e., 6- 
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dependence) of the signal becomes very small, and a model-independent separation of the electromagnetic 
and muon components is impossible. The method of determining the energy scale using the constant intensity 
method is then not applicable. However, the attenuation curve can still be measured experimentally and 
compared with the prediction from simulations, which depends on the energy spectrum of muons produced 
in EAS. In addition, hybrid events with an independent energy measurement allow for a measurement of 
the absolute muon signal normalization with respect to simulations. 

The ratio of muon to electromagnetic signal is determined by three factors in the experimental setup: 
1.) The detector type: thin scintillator detectors have equal response to all minimum ionizing particles and 
thus operate as particle counters. The measured number flux of particles is dominated by electromagnetic 
particles. Shielding can however make scintillator detectors sensitive to muons as well. By contrast, muons 
deposit a large signal in water Chcrcnkov detectors. In this case, the ratio of area to height of the water 
volume determines the EM//i ratio (flatter tanks yielding a larger ratio). 2.) The detector spacing: the 
spacing determines the distance at which the signal is measured [18]. Since the lateral distribution of the 
muon signal is more spread out than the electromagnetic signal, increasing the spacing will increase the 
relative muon contribution in the measured particle flux. 3.) The stage of shower evolution probed: this 
depends on the height above sea level of the experiment and the range of primary energy observed. Showers 
observed very far from the shower maximum have a small electromagnetic component. 

Due to different characteristics, the ground signal will be dominated by the electromagnetic part in some 
experiments (e.g., AG AS A, EASTop, Telescope Array), whereas the muon signal will contribute significantly 
at others (e.g.. Auger, Haverah Park). A possible quantitative criterion for the applicability of the constant 
intensity method is the significance of the signal attenuation observed (e.g., S{6 = 0°)/S{9 = 60°)) with 
respect to the statistical and systematic errors in the signal determination. This criterion corresponds to an 
upper limit on the relative muon signal contribution, which has a very weak dependence on zenith angle. 

8. Discussion and conclusion 

We have shown how Monte Carlo predictions of ground signals can be used to determine the energy scale of 
surface detector experiments, indepedently of the cosmic ray composition and hadronic interaction models. 
This method overcomes the otherwise unavoidable systematics of surface detectors due to the unknown 
cosmic ray composition. In addition, it allows for a clean measurement of the number of muons in extensive 
air showers. In light of the recent detection of a possible GZK feature in the UHECR spectrum, the energy 
scale of cosmic ray experiments is of crucial importance to distinguish between different UHECR source 
scenarios. Hence, it is desirable to determine the energy scale with several methods. The measurement of 
the surface detector energy scale presented here is completely independent of the energy scale determined 
from fluorescence detectors, and contains different systematic uncertainties. 

In this paper, we explored only a single surface detector observable, the signal S{r) at a fixed distance 
from the shower axis. The methodology can be extended to parameterize the signal at difii'erent distances and 
azimuth angles. An extended parameterization like this can then be compared with each detector station 
in a given event, increasing the number of observables for each event. Ideally, perhaps in combination with 
other observables like the rise time [38-40] , this could be used to break the degeneracy of A^^ and energy on 
an event-by-event basis for a surface detector alone. 

It is important to note, however, that air shower universality, the basis of this methodology, can be violated 
by new mechanisms in hadronic interactions in EAS. Recently, the hadronic interaction model EPOS has 
been introduced [41,42]. While the predictions for the depth of shower maximum are within the range of 
the previous models considered here, EPOS shows considerable deviations in the ground signal predictions: 
at 10^^ eV, the EPOS electromagnetic signal seems to be ~ 20% larger than in the other models, while the 
predicted muon signal is 50-70% higher. These differences are due to the production of secondary baryon- 
antibaryon pairs in the GeV range which is strongly increased in EPOS. These baryons then produce more 
muons and a flatter lateral distribution of the signal compared to the other models. These predictions, while 
violating air shower universality, can be constrained by observations using the methodology presented here: 
by separately parametrizing the electromagnetic and muon signals predicted by EPOS, one can infer the 
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relative muon normalization with respect to EPOS which is required by the data. We would like to point out 
that EPOS can be compared with cosmic ray data at lower energies (e.g. KASCADE [43]), as it was done 
in the past with the QGSJetll and Sibyll2.1 models. In addition, accelerator experiments [44] are underway 
to measure the baryon pair production at the relevant energies. One might hope that, once the magnitude 
of the baryon-antibaryon production is understood, hadronic models will converge to a universal prediction 
of the electromagnetic part as shown here for QGSJetll and Sibyll. 

Very generally, the methodology presented here allows for a clean comparison of Monte Carlo simulations 
with air shower data, by separating shower evolution effects from primary composition and high-energy 
interactions. In applying air shower universality, current and future experiments have the potential to tightly 
constrain high energy hadronic models, as well as the energy scale and mass composition of the cosmic ray 
beam. 
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Appendix A. Geometrical effects, ground plane and shower plane signals 

Fig. A.l shows a sketch of an inclined shower hitting the ground. Two detectors at the same distance 
from the shower core correspond to different stages in the shower development, i.e. different DX (Eq. (1)). 
The signal size in inclined showers shows a modulation with ( angle. This modulation or signal asymmetry 
is produced by a convolution of effects [10]. The first one is due to the C dependence of DX for inclined 
showers (Eq. (1)). In addition, there is a geometrical effect which depends on the detector geometry and 
zenith angle of the shower. 

Consider the flux of shower particles $ at a given DX and r, so that the number of particles entering 
the detector is $ • Ad, where Ad is the vector associated with the detector surface with |Ad| equal to its 
area. $ is invariant under rotations around the shower axis, and it only depends on r and DX. For vertical 
showers, the number of particles $ • Aj is invariant under rotations around the shower axis. For inclined 
showers, however, this is not the case, since $ is not parallel to the shower axis (Fig. A.l). We therefore 
expect detectors at the same r but different ^ angles to be hit by different numbers of particles, even if $ 
was independent of DX. 

Hence, the ground-plane signals depend on DX, r and 9. To suppress the dependence on 9 and decouple 
shower development effects from geometrical effects, we define the shower plane signal as the signal generated 
by shower particles passing through the top surface of a detector placed perpendicular to the shower axis. 
This allows us to combine different zenith and ( angles when plotting the shower-plane signal vs DX (e.g., 
Fig. 1). 

While the ground signals can be straightforwardly calculated from the simulation output particles (taking 
into account the statistical weight W due to thinning), the shower plane signals are calculated in the following 
way. 

Consider a particle with weight W and unit momentum vector p hitting a sampling region at ground. We 
will associate a flux $i to this particle in such a way that: 

N^P = W = ^i-Af^ =^iAfPp-S (A.l) 
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Fig. A.l. Sketch of an inclined shower hitting the ground level. The DX of each tank is its distance to the shower maximum, 
projected on the shower axis {DXq it l\ux in the sketch). 



where N'^^ is the number of particles hitting a ground samphng area A'^^ 10^ m'^), and g is the unit 
normal to the ground. Therefore, the flux associated with a particle of weight W is given by: 

W 

P • g AGP 

The number of particles crossing the corresponding sampling area Af^ in the shower plane is given by 
the scalar product of this flux and a vector parallel to the shower axis with a normalization equal to Af^ : 

Nr = '^:-A-r = W%^%=wt4%.-^ (A.3) 

p ■ g ^3 p • g 

where a is a unit vector along the shower axis, and we used Af^ /A'^^ = a • g. 
The shower-plane signal of a detector is then given by: 

gSP = ^ ATf P ^ R^E,,) (A.4) 

i ^ 

where Atank is the area of the top surface of a detector, R{Ep^i) is the detector response to a vertically 
incident particle with energy Ep^i, and the sum runs over all weighted particles falling into the sampling 
area. 
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